Not applicable.
Not applicable.
Not applicable.
Hydraulic fracture mechanics, by far the most popular well stimulation technique, is often plagued by the uncertainties in field parameters for accurate field implementations. For vertical wells, uncertainties in reservoir parameters, such as far-field stresses, may only affect the size of fractures and do not pose many problems otherwise with respect to the geometry of the resulting fracture. However, for inclined (or deviated) wells, additional problems are introduced that cause a significant difference in the geometry of the fracture, both in size and shape, from its designed course, even in the near-wellbore region. Hence, all estimates of fracture behavior and post-fracture production should be made with the knowledge of the highly irregular fracture profile. More often than not, this is not done, causing considerable departures between expectations and reality.
The near-well stress concentration is affected by a number of factors, which include the far field stresses, the well deviation from both the vertical and a plane of principal stress, and the well completion configuration. In effect, the fracture initiation and, consequently, the resulting fracture geometry are greatly influenced by this stress concentration. Incomplete knowledge of all of these factors causes problems during execution of hydraulic fracturing, such as elevated fracturing pressures and unintended screenouts, because of tortuosity, which adversely affects the post-treatment well performance with especially severe effects in high-permeability formations. The uncertainty in magnitude and orientation of far-field principal stresses causes many of the unexplained perturbations in near-wellbore fracture profiles.
The far-field stresses, which are caused by overburden and tectonic phenomena, are supplanted by a new set of stresses when a borehole is drilled. This near-wellbore in situ stress field, in the presence of an arbitrarily inclined borehole, is dictated by the equilibrium equations and depends on the far-field stresses. Stress values are directly related to the state of strains through constitutive equations (elastic, plastic, etc.). When a hydraulic fracture is created at a borehole, the fracture initiation point is important to the fracture propagation, which, in turn, depends on the state of stress around the well. As a result, the presence of the fracture in the formation now redistributes the stresses from their original values without the fracture. In principle, if all of the required reservoir data are known, then the exact fracture profile can be predicted. However, in reality, uncertainty frequently is associated with the reservoir parameters, such as the principal stress orientations and, especially, the magnitude of the intermediate stress. An important consequence is that the resulting fracture geometry will not match its design. More important, in high-permeability fracturing, there is a compelling need to align the well, perforations, and the fracture to prevent or reduce very detrimental tortuosity.
For an open-hole completion, the problem has been studied previously and reported in P. Valko and M. J. Economides, xe2x80x9cHydraulic Fracture Mechanics,xe2x80x9d Wiley, West Sussex, 1995. There are predictive models to evaluate both the fracture initiation pressure and the near-well fracture tortuosity, given the far-field stresses and all the angles that can describe the well position and the fracture initiation point. However, when a fracture is introduced into the formation, no closed form analytical solution is available, and numerical models must be relied on to compute the induced stress profile. Typically, finite element models are used predominantly in such solid mechanics applications.
In many cases, hydraulic fracturing may be performed on a completed well having a casing and sheath. The choice of sheath material, such as foamed cement or neat cement, may affect the fracture geometry significantly due to its material properties. Also, the presence of multiple zones may have other influences in the near-well zone, such as on fracture initiation and fracturing pressure. During hydraulic fracturing of a cemented well, for example, internally pressurized wellbores cause the casing to expand, which induces a tensile stress in the surrounding continuous cement sheath. As a result, the fracture initiation is a function of the cement""s tensile strength and the tensile stresses induced within the cement sheath. However, the effect of the far-field stresses should be included in the field, which is almost always asymmetrical in nature. In effect, both tensile and compressive stresses may act on portions of the cement sheath, thereby making some portions more vulnerable to fracture initiation.
As mentioned, finite element models predominate in such applications. However, finite element modeling can become inefficient and cumbersome for many classes of problems, including fracture mechanics. Finite element models are cumbersome when it comes to complex geometry, in terms of their size, reusability with minor changes, and resources required. An alternative approach, the boundary integral equation method (BIEM), was proposed in the 1950""s for fluid flow applications, and applied in the late 1960""s to mechanical analysis. See, for example, C. A. Brebbia, xe2x80x9cThe Boundary Element Method for Engineers,xe2x80x9d Pentech Press, Plymouth, 1978. The boundary element method (BEM) emerged as a more generally applicable technique during the 1970""s, and has been developed substantially in the following years. See, for example, J. Trevelyan, xe2x80x9cBoundary Elements for Engineersxe2x80x94Theory and Applications,xe2x80x9d Computational Mechanics Publications, Southampton and Boston, 1994. Boundary element techniques are far superior to finite element models, due to ease of use, accuracy, flexibility, and computational speed.
The boundary element method is a numerical technique for analyzing the response of engineering structures when subjected to some kind of xe2x80x9cloading.xe2x80x9d The main feature of BEM is that the governing equations are reduced to surface or boundary integrals only, with all volume integrals removed by mathematical manipulation. Because only surface integrals remain, only surface elements are needed to perform the required integration. So, the boundary elements needed for a 3D (three-dimensional) component are quadrilateral or triangular surface elements covering the surface area of the component. Even simpler, the boundary elements for 2D (two-dimensional) and axisymmetric problems are line segments tracing the outline of the component.
The simplicity of a BEM model means that much detail can be included without complicating the modeling process. In particular, cylindrical holes, such as petroleum wells, can be modeled very quickly, where there is no connection between a hole and the outer surface. Boundary elements also allow analysis of problems that would overwhelm finite element models with too many elements. The system matrix for boundary elements is often fully populated (i.e., dense) and non-symmetrical, but is of significantly smaller dimension than a banded finite element global stiffness matrix.
Because boundary elements are simply lines for 2D and axisymmetric problems, there needs to be a convention used for determining which side of an element is the free surface and which side is inside the material. It is most convenient to use the direction of definition of the element connectivity as the indicator of this orientation. Under this convention, as will be appreciated by those skilled in the art, if the direction of all elements in the model were reversed, we would be modeling the entire infinite universe surrounding a void shaped like the boundary element mesh. In petroleum well applications, these boundary elements are very useful since a few elements can model the problem very accurately where several thousand finite elements likely would be necessary.
The boundary elements are located only on the surface of the component, as are the nodes of the elements. This means that the locations at which the boundary element results are found are only on the surface of the component. It is possible to extract the results for any internal point(s) inside the material simply from the solution over the boundary. The results are not just found by extrapolation, but by using an accurate integral equation technique very similar to that used for the solutions over the boundary elements.
Boundary elements also allow us to define models consisting of a set of sub-models, or zones. Zones are boundary element models in their own right, being closed regions bounded by a set of elements. They share a common set of elements with the adjacent zones. These xe2x80x9cinterfacexe2x80x9d elements, which are completely within the material and not on the surface, form the connectivity between the various zones. This zone approach can be employed when a component consists of two or different materials, when components have high aspect ratio, when elements become close together across a narrow gap leading to inaccurate results, or when computational efficiency needs to be improved.
This boundary element method eliminates the necessity for nested iterative algorithms, which are unavoidable when domain integral methods, such as finite difference methods and finite element methods, are used. Using BEM, it is easier to change a model quickly to reflect design changes or to try different design options. The boundary element method is highly accurate, because it makes approximations only on the surface area of the component instead of throughout its entire volume.
The solution to the forward problem using well known calculations determines the induced stress concentration at a point for known internal pressure and far-field conditions, with or without fracture. It is quite useful in avoiding highly undesirable situations a priori or in determining the ideal location of a new hydraulic fracture. For a well, the natural boundary conditions are specified in the form of traction at the far-field boundary and internal pressure at the wellbore. Once these are known, the geometry of the fracture can be modeled in the well using the method shown in P. Valko and M. J. Economides, xe2x80x9cHydraulic Fracture Mechanics,xe2x80x9d Wiley, West Sussex, 1995. A typical conclusion would be that deviated wells are generally far less attractive hydraulic fracture candidates than vertical wells or horizontal wells that follow one of the principal stress directions.
A brief summary of the development of the boundary integral equations for static stress/displacement problems now is presented. The boundary integral equation for elastostatics is derived from Betti""s Reciprocal Theorem, as will be appreciated by those skilled in the art. The BEM is then derived as a discrete form of the boundary integral equation. The reciprocal theorem states that, for any two possible loading conditions that are applied independently to a structure such that it remains in equilibrium, the work done by taking the forces from the first load case and the displacements from the second load case is equal to the work done by the forces from the second load case and the displacements from the first load case. For example, if the two loading conditions are called conditions A and B, we can write:
ForcesAxc3x97DisplacementsB=ForcesBxc3x97DisplacementsA
Now consider an arbitrary body shape made of a certain material and subject to certain boundary conditions (e.g., loads, constraints, etc.), as shown in FIG. 1. The volume of the body is denoted V, and its surface is denoted S. The tractions, displacements, and body forces are denoted as t, u and b, respectively. Also, define a complementary problem in which the same geometry is subjected to a different set of loads, as shown in FIG. 2. In this complementary condition, the variables are the tractions t*, the displacements u* and the body forces b*. Using the reciprocal theorem, the work done by the forces in the real load case (t,u,b) and the displacements from the complementary load case (t*,u*,b*) are equated to the work done by forces in the complementary load case and the displacements from the real load case, or                     ∫        S            ⁢                                    t            *                    ·          u                ⁢                  ⅆ          S                      +                  ∫        V            ⁢                                    b            *                    ·          u                ⁢                  ⅆ          V                      =                    ∫        S            ⁢                                    u            *                    ·          t                ⁢                  ⅆ          S                      +                  ∫        V            ⁢                                    u            *                    ·          b                ⁢                              ⅆ            V                    .                    
If the body forces in the real load case are ignored, the result is                     ∫        S            ⁢                                    t            *                    ·          u                ⁢                  ⅆ          S                      +                  ∫        V            ⁢                                    b            *                    ·          u                ⁢                  ⅆ          V                      =            ∫      S        ⁢                            u          *                ·        t            ⁢                        ⅆ          S                .            
It is helpful for the complementary load case to represent a type of point force. The form of the point force is the fictitious Dirac delta function. This condition gives rise to boundary reactions, where the component is restrained in the complementary condition, and also a displacement field to consider for the complementary case. The Dirac delta function is defined for all field points y and source point p in the volume V as       Δ    ⁡          (              p        ,        y            )        =      {                                                      0                                                      y                ≠                p                                                                        ∞                                                      y                =                p                                                    ⁢                  
                ⁢                              ∫            V                    ⁢                                    Δ              ⁡                              (                                  p                  ,                  y                                )                                      ⁢                          xe2x80x83                        ⁢                          ⅆ              V                                          =      1      
Because the integral of the Dirac delta function is 1 over the volume V, the volume integral of the Dirac delta function and the real load displacement can be reduced such that             ∫      V        ⁢                  Δ        ⁡                  (                      p            ,            y                    )                    ⁢      u      ⁢              ⅆ        V              =            u      ⁡              (        p        )              .  
Thus, the choice of the Dirac delta function is useful to eliminate the volume integral term in the reciprocal equation. Also, the traction and displacement fields can be estimated (from classical theory) when a point force of this type is applied at a point source p. These are known functions, called xe2x80x9cfundamental equations.xe2x80x9d For 2D problems, the displacement in the complementary load case in the (i, j) direction is given by             u      ij      *        =                  1                  8          ⁢                      πμ            ⁡                          (                              1                -                v                            )                                          ⁡              [                                            (                              3                -                                  4                  ⁢                  v                                            )                        ⁢            ln            ⁢                          1              r                        ⁢                          δ              ij                                +                                    r              i                        ⁢                          r              j                                      ]              ,
where xcexc is the material shear modulus, v is Poisson""s ratio, r is the distance between the source point p and the field point y, and the components of r are ri and rj in the i and j directions. The traction fundamental solutions are given simply as       t    *    =                    ∂                  u          *                            ∂        r              ⁢                            ∂          r                          ∂          n                    .      
Thus, the volume integral term is reduced simply to u(p), and a value of u* and t* for a given source point p can always be calculated, so the reciprocal theorem equation can be rewritten as             u      ⁡              (        p        )              +                  ∫        S            ⁢                                    t            *                    ·          u                ⁢                  ⅆ          S                      =            ∫      S        ⁢                            u          *                ·        t            ⁢                        ⅆ          S                .            
To remove the last non-boundary term in the equation, specify that the point force is somewhere on the boundary and use a constant multiplier c(p)=1 when the fictitious point source is completely inside the material, and c(p)=0 when the point source is on a smooth boundary. Then, the reciprocal equation can be rewritten as                     c        ⁡                  (          p          )                    ⁢              u        ⁡                  (          p          )                      +                  ∫        S            ⁢                                    t            *                    ·          u                ⁢                  ⅆ          S                      =            ∫      S        ⁢                            u          *                ·        t            ⁢                        ⅆ          S                .            
To integrate numerically the functions u* and t*, divide the surface S into many small segments or boundary elements. The integration is then performed over small sections of the boundary surface S, and their contributions are added together to complete the surface integrals. In this discrete form, the surface integral equation may be rewritten as                     c        ⁡                  (          p          )                    ⁢              u        ⁡                  (          p          )                      +                  ∑        elem            ⁢                        ∫          S                ⁢                                            t              *                        ·            u                    ⁢                      ⅆ                          S              elem                                            =            ∫      S        ⁢                            u          *                ·        t            ⁢                        ⅆ                      S            elem                          .            
While the finite order boundary elements, such as constant, linear, or quadratic, etc., are used to provide small areas for numerical integration, the corresponding nodes provide a set of values for interpolation. The discrete form of the boundary integral equation has as its unknowns the displacements and traction distributions around the boundary of the component. This means that when we perform the integrations over every element for any position of the source point, we obtain a simple equation relating all of the nodal values of displacement and traction by a series of coefficients,                               1          2                ⁢                  u          i                    +                                    h            ^                    i1                ⁢                  u          1                    +                                    h            ^                    i2                ⁢                  u          2                    +      ⋯      +                                    h            ^                    in                ⁢                  u          n                      =                                        g            ^                    i1                ⁢                  t          1                    +                                    g            ^                    i2                ⁢                  t          2                    +      ⋯      +                                    g            ^                    in                ⁢                  t          n                      ,
where i represents the ith component of displacement and n represents the number of nodes on the boundary. In doing so, the whole system of equations can be written in the simple matrix form             H      ⁢              xe2x80x83            ⁢      u        =          G      ⁢              xe2x80x83            ⁢      t        ,
where the (nxc3x97n) square matrices H and G are called the influence matrices, and the terms inside them are the influence coefficients. Depending on the boundary conditions specified, the above set of algebraic equations can be rearranged and solved for the remaining unknowns. Having found the values of displacement and traction at the boundary nodes, the solution for the internal points can be calculated using                     u        ⁡                  (          p          )                    +                        ∫          S                ⁢                                            t              *                        ·            u                    ⁢                      xe2x80x83                    ⁢                      ⅆ            S                                =                  ∫        S            ⁢                                    u            *                    ·          t                ⁢                  xe2x80x83                ⁢                  ⅆ          S                      ,
where p is the internal point source.
The calculations at the internal points contain no further approximations beyond those made for the boundary solution. So, as long as an internal point is not so close to the boundary as to make an integral inaccurate, the results there should be just as accurate as the boundary nodal results.
The hydraulic fracturing of arbitrarily inclined wells is made challenging by the far more complicated near-well fracture geometry compared to that of conventional vertical wells. This geometry is important both for hydraulic fracture propagation and the subsequent post-treatment well performance. The effects of well orientation on fracture initiation and fracture tortuosity in the near-wellbore region have been studied and reported in Z. Chen and M. J. Economides, xe2x80x9cFracturing Pressures and Near-Well Fracture Geometry of Arbitrarily Oriented and Horizontal Wells,xe2x80x9d SPE 30531, presented at SPE Annual Technical Conference, Dallas, 1995. These effects indicate an optimum wellbore orientation to avoid undesirable fracture geometry.
As reported in Sathish Sankaran, Wolfgang Deeg, Michael Nikolaou, and Michael J. Economides: xe2x80x9cMeasurements and Inverse Modeling for Far-Field State of Stress Estimation,xe2x80x9d SPE 71647, presented at the 2001 SPE Annual Technical Conference and Exhibition, New Orleans, La., Sep. 30-Oct. 3, 2001, a closed form analytical solution is developed to calculate the stress state within an arbitrary number of hollow, concentric cylinders, with known internal and external pressures. However, far-field stress conditions are assumed to be symmetrical, so that the one-dimensional problem is analytically tractable. The results of the closed form analytical solution now are summarized. Consider n concentric hollow circular cylinders of known internal diameter (ID) ai and outer diameter (OD) bi. These circular cylinders are denoted by indices i, where i=1 refers to the innermost hollow cylinder and i=n refers to the outermost cylinder. Because no void spaces exist between concentric circular cylinders, ai+1=bi. The pressure PO in the innermost cylinder and the pressure Pn outside the outermost cylinder are assumed known. Each cylinder is assumed to behave in a linear elastic manner with known material properties, while the displacement is continuous between cylinders. The stresses "sgr"jk and displacements uj within cylinder i are given by:             σ              rr        i              =                  2        ⁢                  A          i                    +                        B          i                ⁢                  1                      r            2                                          σ              θθ        i              =                  2        ⁢                  A          i                    -                        B          i                ⁢                  1                      r            2                              xe2x80x83"sgr"rxcex8i="sgr"zzi="sgr"rzi="sgr"xcex8zi=0
      u          r      i        =                    α        i            ⁢              A        i            ⁢      r        -                  β        i            ⁢              B        i            ⁢              1        r            xe2x80x83uxcex8i=uzi=0
where xcex1i and xcex2i are functions of each cylinder""s elastic constants:             α      i        =                            (                      1            -                          2              ⁢                              v                i                                              )                ⁢                  (                      1            +                          v              i                                )                            E        i                                β        i            =                        1          +                      v            i                                    E          i                      ,  
vi is Poisson""s ratio, and Ei is Young""s Modulus. The constants Ai and Bi are determined from the pressure applied to the ID and OD of cylinder i:             A      i        =                  1        2            ⁢                                                  b              i              2                        ⁢                          P              i                                -                                    a              i              2                        ⁢                          P                              i                -                1                                                                          b            i            2                    -                      a            i            2                                          B      i        =                                        a            i            2                    ⁢                                    b              i              2                        ⁡                          (                                                P                                      i                    -                    1                                                  -                                  P                  i                                            )                                                            b            i            2                    -                      a            i            2                              .      
The unknown pressures, Pi, between individual cylinders are determined using the requirement of displacement continuity between individual, touching, circular cylinders. Applying the boundary and continuity conditions leads to nxe2x88x921 discrete, linear equations for the nxe2x88x921 unknown contact pressures. With this solution, the stresses and displacements can now be estimated using the constants Ai and Bi.
The above solution works if the far field stresses are known or symmetrical. However, because that is not often the case, it would be helpful if there were a way to quickly and accurately find the far field stresses, the true well departure angle relative to the principal stress orientation, and to use that information to calculate fracture direction geometries in order to find the most useful placement of a hydraulic fracture.
In one aspect, embodiments of the invention feature techniques for determining and validating the result of a fracturing operation by taking advantage of the accuracy and speed of the boundary equation method of mathematics. While on-line pressure monitoring can provide some useful information about the status of a fracturing operation, it is not enough to characterize completely and uniquely the system, and additional information is required, especially for inclined wells. These measurements monitor the fracturing operation continuously and measure the process variables directly, such as well pressure, wellbore surface stresses, and displacements, which can provide useful on-line information to determine the profile of the propagating fracture.
Use of these embodiments also allows designers and users to better select foam cements and other sheathing materials for their projects. Also, using these embodiments to compare the results for a fractured two-zone case against a non-fractured case will help planners to understand the effect of redistributed stress concentration on the well completion.
Embodiments of the invention feature sensors, for example, piezo-electric sensors, to gather data, such as directional stress measurements from a well site, and model the stress distribution in and around the wells, both in the presence and absence of a fracture. If there is a fracture in the formation, the relative location of the fracture can be interpreted by estimating the stress profile before and after a fracture injection test. The embodiments use processes, which, among other abilities, solve inverse elasticity problems. After determining the fracture profile close to the wellbore, selective and oriented perforation configurations can be calculated and performed, which will provide unhindered flow of fluids from the fracture into the well.
In some cases, the effect of far-field stress asymmetry cannot be excluded in the analysis of multiple zone problems, such as in sheathed wells. For this purpose, embodiments of the invention feature the ability to handle such multiple zone systems.